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NUMERICAL SIMULATION USING PARTIAL 
DIFFERENTIAL EQUATIONS FOR MODELING AND 
CONTROL OF MEDICAL RADIATION FIELDS 


Nicolae Marius ROMAN és 
Horatiu Alexandru COLOSI’,; Mircea Valer PUSCA? 


Rezumat. Aceasta lucrare prezinta o noua abordare, bazata pe un model de ecuatii cu 
derivate partiale (EDP), pentru modelarea analogica si simularea numerica a proceselor 
biomedicale cu parametrii distribuiti. Astfel de procese redau penetrarea radiatiilor in 
radioterapie sau in tehnicile de investigare radiologica. Studiul de fata simuleaza 
propagarea prin diferite medii. O caracteristica importanta a metodei propuse este aceea 
de a simula propagarea prin medii neomogene, cum ar fi tesuturile umane. Lucrarea 
ofera un punct de plecare in tehnica controlului automat al dispozitivelor radiologice. 


Abstract. This paper presents a novel approach, based on a partial differential equation 
(PDE) model, for analogical modeling and digital simulation of biomedical processes with 
distributed parameters. Such processes include penetration of radiation during radiation 
therapy or radiological imaging techniques. The present study simulates the propagation of 
radiation through different media. A valuable feature of the proposed method is its ability to 
simulate propagation along inhomogeneous media such as human tissues. It offers a 
starting point for automated control of radiological devices. 


Keywords: partial differential equations, distributed parameters, radiation field, propagation 


1. Introduction 


Medical research and practice offers numerous applications for numerical 
simulation [1-4]. One such application, automated control of radiation field 
intensity generated during radiation therapy, constitutes an efficient method to 
reduce unnecessary radiation exposure for both patient and device operator [5, 6]. 
The aim of our study has been to develop a model for radiation field control at 
various tissue locations investigated or treated using radiological devices. 


The phenomenon of radiological field spreading is represented on Cartesian 
coordinate system (Op; 0q; Or) in figure 1. A possible pattern of spread is 
considered, particularly in relation to these axes and in relation to time (tf). The 
radiation field intensity y(t, p, g, r) is shown in figures 1b, Ic, and 1d, in which: 
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y(t, P. 9.7) = yoolt,s,u)] = Ky - For (0): Fos (s)-u@). (1) 
The following variables and function were considered in a Cartesian space: 


vet p ego er: spate tag tre (2) 


Fos (5) = +1[Fop?(p) + Fog? (a) + For(r)> (3) 


Fig. 1. Representation of radiation field intensity, y(t, p, g, r) in Cartesian coordinate system: 


a. Cartesian system __b. y on Op axis c. y on Og axis; d. y on 07 axis. 
abs (p) abs (p) 
P 4 P. . 
where: Fop(p)=——e ye + eh se YP? (4) 
ew) Pet 
_ abs (q) _ abs (q) 
Foo (q) = 2 é Q am Q2 Os Q2 (5) 
Cpe) OQ) -Q, 
abs (r) abs (r) 
R a R a 
For(r) = : é “ia re Wi 22 | of so) (6) 
Ri — Ry Ry —R, 


The axes (Op) and (0q) define the horizontal plane, in which the field generator is 
located. Field intensity expressed in (1), and axis (Or) correspond to the depth of 
field propagation. The constants in (4), (5), (6) can be approximated by expert 
procedures, namely: 


P, = Pr ; Py =Ap-P,, where: wp =(4+6) andAp>1_ (7) 
Mp(1+ Ap) 
OQ, =—“! __ ; Q)=49-0, where: pig =(4+6) anddg>1 8) 
Hig (1+ Ag) 
r 
R o_o Ry =Apr-R,, where: Up =(4+6)and Ap >1 (9) 


~ ppd Apy 
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where the final abscises (p fo (q fo (r fe correspond to negligible values (for 


example < 0.05). 
For Fop(p), Foo (9). For(r). the abscises (Py ), (dy ), (r;) may show 


perturbations (discontinuities in tissue structure), estimated by: 


abs(s) abs(s) 
S = S a 
Kos (3) = — Se? where: (10) 
5, -S5 So —S, 
s of Syl -S (11) 
ih \ Ar oda Foo | 
Hs -(1+As) 
Ps Exar. 1D AO GES) 
i aekie G7. Wi +0, +R, oa 
+ + 
Hp Ho ER 
Therefore, Fos (s) can also be written: 
_ ES OIS)N 6) =. Le) ee) 
Fos (s) = a, nl a (13) 
l-As 


Dei DoE a 
or: Fas)= Cy: 2°00 +q°+r joe A? +q°+r | (14) 


1 
Lee c, = Hsit4s). Cy aay: C3=T (15) 


where: Co = : 
1-As Sf S 


The results of Fos(s) from (4), (10), (13) and (14) also represent spatial curves 
deformed by multiple degrees of freedom, for example(is), (As) and (s For 


while the “length constants” (P...), (Q...) and (R...) may exhibit “inertia” or 
“attenuation” of radiological propagation in a tissue. 


Function components Fop(p), Fog (q) or For (r) can be approximated by: 


2 
Foyd (V) = Ayg -€ 80d Wa™) (16) 
Ovd vd 


in which the variable (v) can be (p), (q) or (r), by case. 
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The disturbance index (d) refers to the homogeneity disturbance due to tissue 
heterogeneity located on the abscissa(v,), while (A,7) in (16) represents the 
amplitude of heterogeneity illustrated in figure 2.a. 


Fig. 2. Tissue heterogeneity may lead to homogeneity disturbances of field intensity: 
a. Amplitude b. Trajectory in Cartesian space. 


By convenient choice of the parameter (K,,,) a more or less steep (and 
symmetrical) slope can be ensured. In figure 5b, the deformation effect of 
function Fo,q(v) on the right abscissa (vy) over the component ye (v) can be 
seen. If, with respect to time, the spread of radiological field intensity y(tp,q,r) 
was identical on all three axes (Op), (Og) and (Or) in figures 1b, 1c and 1d, then: 
Forp(t) = Forg®) = For, = For © (17) 


which may also be approximated as a customary form of step response, in relation 
to the control signal u(t) applied to the radiological source. Thus: 


(18) 


ty tp 


Fig. 3. Transient phenomenon. 
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The length (t7)of this transient phenomenon (figure 3) reflects a radiation 
propagation inertia, identical on axes (p, q, 7, s) for which ¢>tf,- tends to the 
unitary asymptote. 


For this slow or fast ascending evolution we considered wp =4+6, Ay >1, 
resulting in the following time constants: 


i ir (T, 
T -_ 677 /ak -In| 2 (19) 
My + A,) T,-7, Af, 


1 


where (¢;) is the moment corresponding to the inflexion Fo, (t;). 


It can be seen in figures 1b, Ic, 1d, and 2 that the final values (Pr-O Tp Ve) 


show the depth of the radiological penetration which depends on the control signal 
and on the more or less homogeneous structure of the environment (tissue). 


If this spread inertia is different in relation to the axes, then condition (17) is not 
accomplished, resulting in: 


Forp(t) x Forg (t) 4 For; (t) (20) 
In this case, formally identical to relation (19), each axis will exhibit: 


lof 


T, 5 =——"-; 1 = 177° T, al 
lp We (+ Jip) 2p Tp *\p (21) 
t 
Tg A rT ig (22) 
Lrg A+Arq) 
t 
T; A Top = ArT, (23) 


Mtr ° (l+ Ary) 
where (¢._/), (Ur...) and (Ar...) are specifically established for each axis. 
To transform into Cartesian coordinates it can be shown that: 
ae top thop the Hp Ty +Toq + To, (24) 
T ? T 
Tp t1op +MNq +T2q tir + Tar Tp tlNg +Tr 


finally resulting in: 


(25) 
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which are the equivalent time constants of the Cartesian space, assuming that 
inertial propagations are different along the three axes (p, gq, r). 


Following this method, For(‘) from (1) becomes: 


T a T: ss 
Pye@ 2 es Te 2 ge @ (26) 
T, -T) T, -T, 


which together with F,, (s) from (10), allows the establishment of the radiological 
field intensity from figure (1), 


Yoolt, 5,u(t)] = Ky, + Foot S) Ugo (t), (27) 


where Fi,(t,5) =, (t): Fos(s)and (K, ) is a weighting coefficient. 


Fos(s) | ae a ‘ 


‘3p 


Fig. 4. Fo; and Fos variations. 


For the two functions F,,.(t) and Fos (s) exemplified in Cartesian space in figure 
(4), the relations become obvious: 


t= 0, +t +t, 3 8, =p +a; +r, + 8, =p, tagtn > (28) 
as they were given in figures 1b, lc, 1d and 2. 


2: Analogical modeling of the proposed adjustment scheme 


In figure 5 a block diagram for adjusting the radiological field with components 
modeled by algebraic and differential (29), (30), and partial differential equations 
(31) is presented. 
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—— 
(+ Bq Yoo F410 Yo + Sor Mor + S20 Hag + Mu +42 ° Moa = 
= Ky (Roo so + Boo Ms + Poo ¥sa) 


aa 
Mat >Q >] >{E} fs 


ig 


[hl 


Fig. 5. Block diagram of the proposed adjustment scheme. 
a9 =Wo-M; Km Yoo =M™o+Tn-™| (29) 
Kip-do + K pr: +Kpr- a2 =c, +Tr-C2; Ke-cgo =Ugpt+TRe-uy (30) 
490° Yoo + 410° Yio + 401 * Yo1 + 420° Y20 + 411° Yi1 + 402 * Yor = 
= Ky -(Po.00'“s0 + P1.00°4s1 + P2.00°Us2) (31) 


The indexes in relations (29), (30) and the first parameter from relation (31), 
correspond to the order of derivative with respect to time (f), and the second 
parameter from (31) corresponds to the order of derivative with respect to variable 
(s) defined in relation (2). Considering the example: 


2 2 2 
d*-c _ 9° ~Y00 _ d= Yoo _ dso 
bina: Ws ae 
S 


= 
e dt2 Ot - Os dt 
the radiological field ( y,, ) defined in (26) is converted into electrical signal (m, ) 


by the measurement transducer M from (29). 


The control error (a, ) results from (29), where (wy ) is the reference signal. The R 
controller with PID behavior from (30) is emitting a control signal (c,), applied to 


the electromagnetic field generator E (microwave), also defined in relation (30). 
The radiation emitting unit, E, generates the execution signal (u, ) representing the 


incident radiological field which is applied to the propagation environments (e.g. 
air, tissues) together with a possible resultant perturbation signal (to), resulting in: 

Usy =Ug + to (32) 
In this case (#) corresponds to the homogeneity perturbation Foyq(v) from (16), 
where v = Sg found in (28). 


Thus, the complex phenomenon of radiological propagation is estimated by partial 
differential equation (31). In the left member of this equation we introduced the 
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approximated solution (27), resulting in the right member: 


Dig — Gy hig Py eg Fy Og ag Ge Ee od, (33) 


02 02 
Po1.00 = G10 * Fog + 2° G9 * Fig + 1° Fo (34) 
Pr.99 = 29 * Fog (35) 
in which partial derivatives ( F,..);(F,,) and (F 4 ).are obtained from relation (27). 
Only signal (u,) is included in the block diagram. The disturbance (i) is an 
independent external signal. 


From the equation system (28)...(31), the matrix with partial derivatives 
of the state vector (Mqp,) represented in (36) has been deduced for the entire 
regulation system in figure 5. In (36), the state vector x has the following 


transposed form: 
af 
x =P} {o}1}4ol Yoo} 10) 


Because the signals m = m(t), c = c(t) and u = u(t) are functions only of time, all 
their partial derivatives with respect to (s) are equal to zero. Details concerning 
the preparation and analogical modeling of systems through (u ae may be found 


in [7, 8]. In our case, (Mqpx) from (36) is composed by: 
x(8 x 1) state vector, 
x5 (8x6) state vector derivative with respect to variable (s), 
x7 (20x1) state vector derivative with respect to variable (1), 


XTs (20 x 6), state vector derivative with respect to variables (f) and (s). 


3. Numerical simulation by (M \,<} and Taylor series 


Numerical simulation requires knowledge of initial conditions ClI(fo,s), for 
Xcy = X(tg, 8) =Xz_1, from which results by partial derivative with respect to (s): 


Xs n= Xs (to.S).=XS KA 


in which sequence (k—1) corresponds to the moment f,_; =(k—1)- At, and 
sequence (k) corresponds to t, =k-At, where (Af) is a sufficiently small 
integration step. 

Using (xe7) and (xs cr) above, we computed using specific operations based on 


symbolic derivatives by indices, the composing elements of (ee) and (rs) 
resulting in: 
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xrcr =Xr(to.s)=x7 44 and xps.c7 = X7s(to.8)= x75 K-41. 


is |S |S |o |e |B |B 


het 
B 


— 
Lo 
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After these preliminaries regarding the initial conditions, we approximated 
through iterative calculation steps, 
6 At 6 rt 
Xe = Xp t epg Xe =Xshat > — Arse = ©=©— BD 
T! T! 
= . T=1 
At each integration step (At) we operated with a number of 20 Taylor series, of 
which 6 for (Mox SIM COk > C1k UOk > Uk ) and 14 for ( YOOk > YOlk 2+» YO6K ). 
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It can be observed that 8 elements (signals) included in Taylor series are from the 
state vector component (x) and the other 12 are from the matrix component (x S ), 


each one included in (M ee). 


7 
: : : : ; At 
The truncation error at each integration step (At) is proportional to (=) , for 


At® 
(M,C,U, You), aS Well as to aa , for ( V9): 


4. Case studies 
4.1. Preliminaries 


a) The radiological field spread has been estimated by the following parameters: 
By Op axis: 
t of = 10; LT = 41; At p = 1,5; Tip = 0,975; Ty» = 1,463 
jee Pa = 9,3 Nidajes Aga = 50; K nq = 10; Hp =4 
Ap =1,1; P, =0,595; P, =0,655 


By Og axis: 
tof = 22 Eg =43: Arq ley Tg = 1,033; ‘by, =lLi7 


qa da = 0B aa;z; Agat=45; Kgq =12; Ho =4,2 
Ag A,4; Q, =0,794; Q, =1111 


By Or axis: 

ty = 15; Upp =45; Apr =L9; T,, =1,149; T5, = 2,184 

rf = 50; Yd 0,37" FF; A,q= 955; K,q = 15; Lr = 4,5 

Ar =1,5; R, =4444 Ro = 6,666 
According to (24, 25), it resulted: 

H, = 4,321; A, =1,711; FT, =1,85 and T, =3,16. 
According to (11, 12), it resulted: 
Ms =4.416;A5 =1,445; S$, =4,71;S, =6,81. 
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b) The structure parameters, and the equation coefficients of the system (29, 30) 
were: 


Ky =|: Ty =0.055; Kp =K, =1; Tr =T, =0,06; Tp = 0,05 


2th, % ! T -Tp 
2-Kyry-Ts 


ee ey ge ee 
TR Kp, Ts DR 2K, Ts (38) 


K pr 
By convenient choice of ( Kz, ), a flexible controller PID for relation (29) shown 
in Tables 2 and 3 has been ensured. 


The intensity field X-ray reference signal ( w, ) and the perturbation signal (tf) in 


figure 5, present a continuous component on top of which periodic components 
may overlap. 


c) The analogical model for the spread of the radiological field expressed by the 
second order partial differential equation in (31) presents the coefficients: 


ago = L a0 = T, +T5 = 5:01 a0 =T, -T> = 5,848; 
ao] = 8; + $5 = [le 
ag2 = S$, So = 32,08; a, =; +T>)-(S; + $5) ee 


The proportionality coefficient K, =1, may be altered over large limits, since the 


variables presented in figures 1b, Ic, Id, 2, 3 and 4, are asymptotic or non- 
asymptotic between 0 and 1. 


4.2. Case study 1: System behavior in open loop 

The control system in figure 5 may be considered in stopped reaction, 
respectively, m,=0. 

The R controller and the radiological field generator ensure a unitary transfer with 
a negligible delay of the reference signal wp = dp =co = Ug. 

The radiological field E generated by the control signal (u,) develops inside a 


propagation space an intensity yoo[t,5,Us9(t)], representing the solution of partial 
differential equation (31). 


Table 1 presents the evolution of the intensity of the above radiological field 
(Yoo ), With respect to time (t), with w, =10, for 11 points from a Cartesian space, 


respectively (0,0,0),...( Dysdpstp)- 
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Table 1. Intensities of the simulated radiological field in open loop 


Dp q r t= 10? 1 3 5 7 9 12 15 18 22 

0 0 0 10° | 0,65 3.47 5.98 7.68 8.71 9.48 9.79 9.90 | 9.965 
O.1-pr | O.1-qr 0.1 -tf 10° | 0.51 2.68 4.64 5.95 6.75 7.34 7.59 7.68 TAS 
0.2-pr | 0.2-q¢ 0.2-t, 10° | 0.43 2.29 3.96 5.08 5.76 6.27 6.48 6.56 6.60 
0.3-pr | 0.3-q¢ 0.3 -tf 10* 1.63 8.64 14.9.) 19.14 21.7 23.6 24.4 24.7 | 24.84 


0.4-pr | 0.4-q¢ 0.4-1, 5:10° 0.123. | 0.653 | 1.127 1.44 1.64 1.78 | 1.845 1.87 1.88 


O.5-pr | 0.5-qr 0.5-1, 107 | 0.044 | 0.233 | 0.403 | 0.517 | 0.586 | 0.638 | 0.660 | 0.667 | 0.671 


0.6-pr | 0.6-q¢ 0.6-r¢ | 410° | 0.022 | 0.115 | 0.198 | 0.255 | 0.289 | 0.315 | 0.325 | 0.329 | 0.331 


0.7-pr | 0.7-q¢ 0.7%, | 2:10% | 0.010 | 0.056 | 0.097 | 0.124 | 0.141 | 0.153 | 0.159 | 0.160 | 0.161 


0.8:p; | 0.8-g¢ 0.8-r¢ 10% | 0.005 | 0.027 | 0.047 0.06 | 0.068 | 0.074 | 0.076 | 0.077 | 0.078 


0.9-p; | 0.9-g¢ 0.9-r¢ | 510° | 0.002 | 0.013 | 0.022 | 0.029 | 0.032 | 0.035 | 0.036 | 0.037 | 0.037 


Pr qr tr | 2:10° | 110° | 610° | 1-107 | 0.013 | 0.015 | 0.017 | 0.0176 | 0.0178 | 0.186 


K, =0,K =1,;w, =10;K,, =1,K,, =0;K,, =0;K,. = De = 2) 
It may be observed that in the near vicinity of the radiation source application 


point, (0,0,0), the radiological field intensity (y),) has a maximum. value, 


presenting an exponential evolution which tends asymptotically towards the 
references value wo =10. The further away we move from the source, towards 


the final penetration depth (py,qf.ry), the more the radiological field intensity 


( yoq ) decreases. 


The (0,3-p¢;0,3-¢ 7 30,3-rg) point corresponds to a structure perturbation 


(Pa>a>'a)» for which the intensity (yoo) presents a remarkable discontinuity, in 


conformity with figure 2, relation (16) and figure 4. It may arbitrarily be 
considered that this discontinuity presents a “zoom effect’, by concentrating 
radiological field intensity (yoo). Of course the opposite may also be considered, if 
in (24) the amplitude A, <0. 


4.3. Case study 2: System behavior in control loop 


In this case, the transducer M of radiological field intensity (v9 ) achieves a 


control error a= w—™m, using the measured signal (m) and a reference signal (w) : 


which is processed by the R controller, with PID structure. Through the transfer 
function of the regulator, 
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1 
K,j+—-K,+sK 
Le PR 5 IR DR 
= (39) 
La Lars 


the three effects: proportional K pr, integration K7p and derivative K pr, can be 
modified over large limits if the K ;y coefficient is wisely chosen from (38). 


Table 2 presents the evolution of the radiological field intensity (yoo), with 
respect to time (t) , with w=5, for the same 11 points of the Cartesian space. By 
convenient adjustment of the controller R, an intensity yog =5 was ensured in 
steady state regime. 

For Kgg=15)K,=1; wSa3 tf =22; Kry conveniently chosen to 
ensure yo9 =5 in steady state regime, the evolution is non-periodically ascending 
towards yog =5, at tr = 22. The values of curves are identical, which means 
that the numerical simulation stages have been performed correctly. 


The only exception corresponds to the perturbation _ point 
(0,3 Py 0,3q Fs 0,3 lf ), where the “zoom effect” of ry ) presents a slight peak, 


at 6.45, then returning to the same asymptote yoo =5. 


In this perturbation regime, the relatively small value of Ky =0,282 constrained 
an important growth of the PID controller effects. 


Table 2. Intensities of the simulated radiological field in control loop 


Pp q r yoo | Kex | t=107 | “1 3 5 i 9 12 15 18 23: 


0 0 0 Yoos | 15 6-107 | 0.33 | 1.76 | 3.03 | 3.88 | 4.38 | 4.76 | 4.91 | 4.97 | 4.999 


O.L pp | OL ag, | OL 7% | yoos | 11.6 6-107 | 0.33 | 1.76 | 3.03 | 3.88 | 4.39 | 4.76 | 4.91 | 4.97 | 4.999 


0.2 pp | 0.2q, | 0.217 | yoo | 11.4 | 1-10 | 0.33 | 1.756 | 3.02 | 3.87 | 4.39 | 4.77 | 4.92 | 4.98 | 5.002 


0.3 py 0.3 a | 0.3 r% | yooo | 0.282 | 2-104 | 2.28 | 5.62 6.45 | 6.33 | 5.99 | 5.54 | 5.27 | 5.13 | 5.05 


0.4 py 0.4 aq, | 0.477 | yoor | 2.25 5-10° | 0.33 | 1.758 | 3.03 | 3.88 | 4.38 | 4.76 | 4.91 | 4.97 | 4.998 


0.5 py 0.5 a, | O57 | yoor | 1.01 6-10’ | 0.33 | 1.76 3.03 | 3.88 | 4.38 | 4.76 | 4.91 | 4.97 | 4.997 


0.6 py 0.6 q | 0.6 r¢ | yoog | 0.500 | 6-107 | 0.33 | 1.757 | 3.02 | 3.87 | 4.38 | 4.76 | 4.91 | 4.97 | 4.996 


0.7 pr 0.7 gy | 0.7 r¢ | yoou | 0.245 | 6-107 | 0.33 | 1.751 | 3.02 | 3.87 | 4.37 | 4.76 | 4.91 | 4.97 | 4.999 


0.8 pr 0.8 gq | 0.8 r¢ | yoor | 0.117 | 6-107 | 0.33 | 1.765 | 3.03 | 3.88 | 4.39 | 4.77 | 4.92 | 4.97 | 4.999 


0.9 py 0.9 qr | 0.9 r¢ | yoor | 0.057 | 6-107 | 0.33 | 1.75 3.01 | 3.86 | 4.37 | 4.75 | 4.90 | 4.97 | 4.994 


Pr Ch If yoox | 0.027 | 6-107 | 0.33 | 1.76 3.03 | 3.88 | 4.39 | 4.76 | 4.92 | 4.97 | 4.998 
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Thus, a constant value of field intensity at yoole f = 5 in steady state regime, for 
reference w=5, over all 11 points yog4, Yoog>---» Yoox has been ensured. 


Table 3. Influence of the Kgy coefficient on the R controller 


00 K EX KPR KIR Kpr voles) us(ty | 
Yooa 15 1.0125 0.202 1.181 4.999 5.009 
Yoon 11.6 1.309 0.261 1.527 4.999 6.464 
yooc 114 1.332 0.266 1.554 5.002 10.68 
Yoop 0.282 53.86 10.74 62.84 5.05 19.69 
YOOE 2.25 6.750 1.345 7.876 4.998 37.51 
Yoor 1.01 15.037 3.00 17.546 4.997 74.45 
Yoo 0.500 30.37 6.06 35.44 4.996 150.71 
YooH 0.245 62 12.37 72.33 4.995 308.94 
Yoor 0.117 129.8 25.9 151.46 4.999 639.7 
Yooy 0.057 266.4 53.16 310.9 4.994 1329.8 
Yoo 0.027 562.5 112.2 656.35 4.998 2781 


The further the M transducer is placed from the radiological field-generating 
device, the more intense the three effects of the controller, resulting in larger 


K,,»K, and K,,, coefficients. 


PR? 


Near the structural discontinuity, corresponding to the yogp concentration, the 
controller ensures the attenuation of the resulting perturbation, by constraining 
larger values for the three PID effects. 

In a stationary regime yoole f 2 5, the execution signals us (r f) present roughly 


uniform increasing values, for progressively increasing distance of the transducer 
M from the field generating device. 
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Conclusions 


The paper studies a modeling approach of a radiological propagation process, 
based on a system of three differential equations and a partial derivative equation 
(PDE) with three space variables and one time variable. 


For analogical modeling it operates using the matrix of partial derivatives for the 
state vector M 4,,, associated with Taylor series for computational simulation. 


A controller (R) can be tuned in order to ensure the control of radiological field 
intensity at different points where field sensors and transducers are located. 
Multiple degrees of freedom can be ensured using structural parameters of the 
model and external elements such as the above mentioned controller, transducers 
and sensors. Altering of environmental propagation can be distinguished using 
structural parameters (ue. ) and Gi ), in time as well as in space. 


The model accounts for the presence of structural tissue discontinuity. Specific 
time-space dynamics are taken into account for each space axis (O,,,0,,0,). 


Units of measurement for length (p.4, r), time (t), execution signals (u), 
radiation intensity (y) and other parameters may be chosen as needed, but will 
only be valid for each individual case study. 


The PDE model may be adapted for various processes of radiation or thermal 
propagation, with applications in both medical diagnosis and treatment. 
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